Random Field Potts model with dipolar-like interactions: hysteresis, avalanches and 

microstructure 
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A model for the study of hysteresis and avalanches in a first-order phase transition from a single 
variant phase to a multivariant phase is presented. The model is based on a modification of the 
Random Field Potts model with metastable dynamics by adding a dipolar interaction term truncated 
at nearest neighbors. We focus our study on hysteresis loop properties, on the three-dimensional 
(3D) microstructure formation and on avalanche statistics. 
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I. INTRODUCTION 

Microstructure formation in first-order phase transi- 
tions is a phenomenon that has been studied by physi- 
cists, mathematicians and engineers B It is im- 
portant not only from a fundamental point of view, but 
also for applications, due to its relationship with material 
properties. Microstructures occur in ferroic systems (fer- 
romagnetic, ferroelectric and/or ferroelastic) which are 
driven through a first-order phase transition (FOPT) in 
which some symmetry operations of the parent phase are 
lost. The product phase (usually less symmetric) may 
appear in energetically equivalent variants which are re- 
lated by the symmetry operations that have been lost in 
the transition. The obtained microstructures correspond 
to the arrangement of such equivalent variants and are 
decided by the interplay of many energetic terms: inter- 
face energy, surface energy and long-range interactions. 

Until now many of the studies on microstructures have 
focused on the determination of the optimal variant con- 
figuration minimizing a certain thermodynamic potential 
that takes into account the above factors and external 
conditions. Nevertheless, in many cases, when real ma- 
terials are studied, such optimal microstructures are not 
observed. This is mainly due to two important factors: 
(i) the existence of disorder sources of very different na- 
tures both in the bulk and on the surfaces and also (ii) 
the athermal character of the phase transition dynam- 
ics: when the temperature is not very high, the ener- 
getic barriers that separate the optimal solutions from 
the parent phase cannot be overcome. Thus, the system 
evolves following metastable paths which locally optimize 
the system energy, but are far from the trajectories ob- 
tained from a global minimization principle. An interest- 
ing suggestion on the behavior of microstructure forma- 
tion comes from the glass-jamming transition framework 
(see Refs. 0, H and references therein), which is asso- 
ciated with the so-called kinetically constrained model. 
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These models are stochastic lattice gases with hard core 
exclusion, with the addition of some local constraints, 
which mimic the geometric constraints on the possible 
rearrangements in physical systems. Similar behavior is 
quite likely to arise in microstructure formation. 

The use of continuum models derived from elasticity 
theory [1, 0, B H, has been proposed as another ap- 
proach to microstructure formation. Some of these mod- 
els have been successful in explaining microstructures, 
hysteresis and avalanches. Nevertheless, they are very 
time consuming from a computational point of view. For 
this reason, in many cases only 2D problems have been 
addressed, and even this item presents difficulties con- 
nected with large statistics and with the scanning of the 
model parameters in order to study their influence. Our 
aim here is to find a statistical mechanical lattice model, 
easy to simulate and which allows for the study of statis- 
tics of microstructure sequences dynamically generated 
in athermal systems, under the influence of disorder. 

The Random-Field Ising model (RFIM) [ill with 
metastable dynamics is one of the simplest models for 
the study of the combined effects of disorder and ather- 
mal evolution. It is formulated in a magnetic language 
for a spin reversal transition, driven by an external mag- 
netic field H . The only degrees of freedom are spin vari- 
ables defined on a lattice [i = 1, N), which take values 
Si = ±\ on the i-ih lattice site. The RFIM enables com- 
putation of hysteresis loops M{H) corresponding to the 
behavior of the order parameter M — St as a func- 
tion of H as well as the analysis of the intermediate states 
between the negatively saturated initial state {Si = — 1} 
and the final positively saturated state {Si — +1}, and 
vice versa. In particular, the RFIM has been successful 
in understanding the avalanche dynamics (Barkhausen 
noise [l^) that joins the intermediate metastable states 
and shows absence of characteristic scales for a critical 
amount of disorder. Moreover, the RFIM displays sev- 
eral other interesting properties [ill, HSl' it exhibits a 
well-defined rate-independent trajectory, it shows return 
point memory, it satisfies the abelian property and, from 
a computational point of view, it is fast to simulate tra- 
jectories in relatively large systems [l3 |. 
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Nevertheless, the usefulness of the RFIM for the study 
of microstructures is almost null. This happens because 
the parent and product phases in the RFIM are the pos- 
itively magnetized phase and the negatively magnetized 
phase. These two phases are single variant and totally 
equivalent from a symmetry point of view. Consequently, 
the obtained hysteresis loops are symmetric under the 
exchange H —H and M —M, there is absence of 
latent heat associated with the FOPT and the domains 
are spherically symmetric (except for some short-range 
correlations due to lattice symmetries). 

Within this framework, the aim of the present work 
is to explore some modifications that should be intro- 
duced in the RFIM in order to obtain 3D microstruc- 
tures without losing, as far as possible, some of the use- 
ful RFIM properties that we have mentioned above. A 
first step in this direction was done several years ago 
by defining the Random Field Blume-Emery-Griffiths 
model, in which the 'spin' variables take three different 
values {Si = —1,0, 1) with metastable dynamics In 
this case, the FOPT takes place from a single variant par- 
ent phase, represented by Si — 0, and a product phase 
with two variants Si — +1 and Si = —1. In the cited 
work, the hysteresis loops, phase diagram and avalanche 
distribution were studied for this type of simple case. In 
the present paper we would like to go one step forward. 
This will be done by starting from the Random Field 
Potts model [l^ with metastable dynamics. In the Potts 
model the spin variables can take an arbitrary number of 
values. This model allows phase transitions to be stud- 
ied from a non-degenerate phase Si = to a multivariant 
product phase. We will explore the effect of an extra in- 
teraction term (of a dipolar nature, but truncated to the 
nearest- neighbor approximation), which will be necessary 
in order to produce microstructures (for the introduction 
of dipolar interactions in RFIM, see [l3|)- It is not our 
aim to focus on the detailed analysis of any particular 
transition. We will study a model that, from the point of 
view of symmetry, would correspond to a transition from 
a cubic phase (single variant) to a tetragonal phase (with 
three equivalent variants), although neither the detailed 
interactions nor the external constraints will be tuned for 
the particular modeling of such transitions in ferroelastic 
systems. This will be the aim of a future work [2^ . 

The paper is organized as follows: in section |TT] we 
introduce the hamiltonian of the model. In section IIIII 
we detail the metastable dynamics that has been used 
for the simulations. Section |W] is devoted to the discus- 
sion of the obtained results: in section HV Al we present 
our analysis on the shape of hysteresis loop cycles as a 
function of the model parameter values. In this section 
it is shown that the loops happen to be unsymmetric, 
in contrast to the RFIM results. The microstructures 
are analyzed in section ITVBl where we discuss three dif- 
ferent regimes corresponding to different ranges of the 
parameters values. Moreover, in section flV CI we present 
the statistical analysis of avalanche behavior. Finally, we 
summarize and discuss the future perspectives in section 
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II. MODEL 

The model can be defined on any regular lattice. We 
will consider a simple cubic lattice of size N = L x L x 
L with periodic boundary conditions. At each lattice 
site we define a variable Si {i — 1,...N), which can 
take four different values that we will call 0, x, y and z. 
We can choose different representations for our variables 
but it is convenient to consider a vector Si having three 
components: we will indicate the four possible values as 
= (0, 0, 0), a: = (1, 0, 0), y = (0, 1, 0) and z = (0, 0, 1). 

The order parameter for the phase transition under 
study here is M = '^i{Si)'^, where the sum spans over 
the whole lattice. M represents the amount of the sys- 
tem that has transformed from the cubic to the tetrag- 
onal phase. By following the analogy with the magnetic 
case, we will refer to M as the total magnetization of the 
system. Moreover, we define the normalized magnetiza- 
tion m, as m = M/N. We will drive the system by an 
external field H coupled to AI, since we are interested in 
the transition from the phase to the multivariant phase 
which will be composed of regions (variants) in the states 
X, y and z. The field H would correspond to the driving 
effect of the temperature in athermal structural transi- 
tions. We will start by decreasing i7, from the M = 
state. We consider the following hamiltonian: 

NN NN ,^ \(Q ^ ^ 

n = -kY, 60., s,) + A ^ • '^'if ' 

<ij> <ij> ' 

N 

^HY,{S,f + Hd.s (1) 

i 

The first term is a Potts exchange term extending to 
nearest-neighbor (n.n.) pairs. The parameter k will be 
always positive, favoring the ferromagnetic interaction 
between spins in the same state. In the following, without 
loss of generality, we will consider k ^ 1. 

The second term is a dipolar interaction truncated to 
n.n. pairs. As we have already said, the aim of adding 
this term is to generate some simple microstructures. To 
include higher order terms would lead to more realistic 
ones. The vector Tij is the lattice vector joining 'spins' 

Si and Sj . We will study the cases with A < and A > 
separately. As can be easily seen from Eq. [1] in fact, the 
A < case corresponds to favoring the growth of prolate 
(needle-like) domains parallel to the spin direction. On 
the other hand, for A > such a growth is not favored, 
but as it is partially compensated by the exchange term 
it basically corresponds to the formation of oblate (disk- 
like) domains, perpendicular to the spin direction. We 
will illustrate these features in Sec. IIVBI 

The third term of the hamiltonian accounts for the 
interaction between the system and the external field H . 
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This field will be driven from very high positive values 
to very negative values (and vice versa) in steps Ai?, 
mimicking an adiabatic triangular driving force (i.e. field 
frequency a; — > 0) . By a deliberate abuse of language we 
will refer to the step Ai? as the driving rate. One can 
notice that it is possible to add a second driving term 
Si which would be convenient for the study of the 
transitions from one variant to another, mimicking, for 
instance, the effect of an applied external stress. 

The last term Tidis accounts for the quenched disorder 
of the system. We will restrict ourselves to a random 
field type with zero averages. However, there are still 
several possibilities for such a hamiltonian term because 
the random fields can couple either to Si or to the order 
parameter (S'i)^. We will thus consider: 
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where gi is a three-component vector random field, whose 
components are extracted from a gaussian distribution 
A'^(0, 1) with zero mean and unitary standard deviation, 
and hi is a scalar field, again extracted from a gaussian 
distribution A^(0, 1). The parameters a and p control the 
amount of quenched disorder in the system. 

In order to compare our model with the standard 
RFIM, we define the total amount of disorder o-g = 
+ p^. In practice, the two disorder terms can be un- 
derstood as arising from a local random field fi, whose 
components are correlated, being 

/i = (^{gix,giy,g%z) + p{h.,, hi, hi), (3) 

so that {fl) = Ul)^ = ifl) = and = 

ifixfiz) — ifiyfiz) — P ■ 



III. DYNAMICS 

There are several possibilities for the choice of 
metastable dynamics. In Fig. [T] we show examples of 
hysteresis loops obtained with three possible choices of 
dynamics. At first sight, the three loops look very simi- 
lar. In all the cases we start from a metastable state, we 
increase or decrease the field by a AH step and then, at 
constant field, we recursively minimize the system energy 
by using a local rule based on single-spin changes. Only 
after a new metastable state is reached we proceed with 
a new field change AH. 

In the extremal selection + extremal update case, we 
scan the whole system and check which variable Si can 
change to a new value with the minimum energy differ- 
ence ATi.. The proposed change is accepted if this mini- 
mum ATi is negative. In the random selection + random 
update case, we randomly choose a spin on the lattice and 
propose a random change to a new value. If the proposed 
change implies ATI < the change is accepted. Finally, 
in the random selection + extremal update dynamics, we 



FIG. 1: (Color on line) Hysteresis cycles for three different 
dynamics as explained in the text. The parameters of the 
simulations are: Aff = 0.05, L = 16, A = — 10, a = 5 and 
p = 2. In the inset: a magnification of a loop region: the 
magnetization values for the three dynamics coincide only for 
some field values. 



randomly choose a spin on the lattice and check among 
the three possible new values which one represents a min- 
imum AH.. If this minimum value is negative we accept 
the change. From a computational point of view the first 
choice is much more time consuming than the other two, 
since the effort scales with . 

Although the loops are very similar, detailed analysis 
reveals that the obtained hysteresis loops, as well as the 
sequence of metastable states, are not identical. This 
tells us that the proposed model is not abelian and that 
the final state will depend on the order in which unstable 
spins will be changed. In order to ensure some robustness 
of the results, we are thus forced to choose extremal selec- 
tion + extremal update dynamics, that is: to propose the 
optimal spin change among the whole lattice and among 
all the three possible final values at each time step. This 
kind of dynamics is deterministic and thus, by definition, 
independent of the updating order. We will keep to this 
dynamics for the rest of the paper. 

We now study the effect of changing the value of AH . 
In Fig. [2] we show three hysteresis cycles obtained for dif- 
ferent values of the driving rate AH using the extremal 
selection+extremal update dynamics. A detailed analysis 
reveals that the differences between the three loops can 
be attributed to the fact that driving with a smaller AH 
allows more metastable intermediate states to be found, 
but for the same applied field values, in the three real- 
izations not only the magnetization, but also the final 
microscopic configurations reached are the same. The 
independence from the field rate is an important prop- 
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FIG. 2: (Color on line) Hysteresis cycles for various val- 
ues of the applied field rate /S.H with the extremal selec- 
tion+ extremal update dynamics. The parameters of the sim- 
ulations are: L = 16, A = —5, a — 4 and p = 0. In the inset: 
a magnification of a portion of the hysteresis cycle. 



erty from the point of view of the simulations, since it 
allows us to use a relatively large AH for the study of 
the properties of hysteresis loops. 




FIG. 3: (Color on line) Examples of hysteresis loops for dif- 
ferent values of the parameter A: (a) A — —1, —5, —10, —15 
and (b) A = 1,5,10. In aU the cases, L = 16, AH = 0.01, 
(T = 3 and p — 



IV. RESULTS 

We have performed numerical simulations of systems 
with sizes L = 8,16,32,40 and 60, averaging over 
10^ — 10'^ realizations of the quenched random fields. We 
have focused our analysis on hysteresis loops behavior, 
on microstructure formation and on the statistical prop- 
erties of the avalanches. 



A. Hysteresis loops 

In Figs. [21 |3]and[5l we show some examples of hystere- 
sis loops simulated with our model in order to illustrate 
the effect of the different hamiltonian parameters. 

We consider the cases with A < (Fig.[3ja)) and A > 
separately (Fig. El^b)), because they show a clearly dif- 
ferent behavior. In the first case, which corresponds to 
the formation of prolate domains, the more negative A 
is, the larger the width of the loop. In the case of very 
negative values, the loops start to exhibit a plateau in 
the increasing field branch: the retransformation to the 
phase is done in two separate steps. (This effect will be 
discussed below). For the second case, increasing lambda 
towards positive values increases the tilt of the hysteresis 
loop so that saturation in the transformed phase can only 
be obtained when the field is very negative. This effect is 



due to the competition between the Potts and the dipolar 
terms. Many domains in the final stages of the transfor- 
mation are frustrated and can only be transformed by a 
very negative , as occurs in ferromagnets that contain 
a small percentage of antiferromagnetic bonds. 

In Figs, m and [5] we show the effect of the two dis- 
order parameters a (Fig. 21 in the low (a) and high (b) 
p regimes) and p (Fig. [5l in the low (a) and high (b) a 
regimes). In all the cases it can be seen that increasing 
the amount of disorder increases the tilt and decreases 
the width of the loop. Moreover, as expected, for low val- 
ues of the amount of disorder (cr or p) the loops exhibit 
sharp discontinuous (ferromagnetic-like) behavior. This 
feature is in agreement with recent results on the stan- 
dard RFIM, concerning the observation that the transi- 
tion from sharp to smooth loops can be induced by dif- 
ferent kinds of disorder parameters: not only the random 
field variance a but also random anisotropy flS*! , the va- 
cancy concentration [19], etc.. Our model shows that the 
correlation with the random fields of intensity p can also 
act in a similar way. 

Let us now discuss the plateau observed in Fig. ^a.) 
in the increasing field branch. As shown in the examples 
in Fig. [6l this plateau occurs at smaller magnetizations 
when the system size is increased. This suggest that it 
may be due to the stabilization of "slab" domains that 
cross the whole system from one face to the other and 
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FIG. 4: (Color on line) Examples of hysteresis loops for var- 
ious values of the disorder parameter a: (a) for a = 4,5, 6, 
with p = 0; and (b) for a = 2,4, 10 with p = 5. In all the 
cases A = -3, L = 16 and AH = 0.01. 



that due to the periodic boundary conditions behave as 
infinitely large. Such slabs become less and less frequent 
by increasing the system size. This suggestion has been 
confirmed by analyzing sequences of configurations. An 
example will be discussed in Sec. IIVBI 

In order to perform a quantitative analysis of the hys- 
teresis loops it is important to measure some of their 
properties. One of the most studied hysteretic features 
is loop area. In fact, it represents the amount of en- 
ergy dissipated during a cycle and thus it is an important 
quantity to be controlled both from the theoretical and 
materials application point of view. In Fig.[7]we show the 
loop area, averaged over several disorder configurations, 
as a function of the parameter A, for various values of a. 

As can be seen, the area shows a much more impor- 
tant dependence on A for negative than for positive A. 
This behavior can also be seen by studying the coerciv- 
ity (amplitude) of the hysteresis cycles at m = 0.5, which 
displays a very similar dependence on A. 

In Figs. [21 m and [5] we can see that the hysteresis cy- 
cles obtained with our model are asymmetric, i.e. the 
decreasing field branch (transformation) cannot be re- 
lated by an inversion operation to the increasing field 
branch (retransformation) . This is an interesting prop- 
erty since, experimentally, materials displaying a transi- 




FIG. 5: (Color on line) Examples of hysteresis cycles for vari- 
ous values of the disorder parameter p: (a) for p = 1, 2, 4 with 
a = 1.5; and (b) for p = 1, 5, 10 with a = 4. In all the cases, 
A = -3, L = 16, and AH = 0.01. 



tion to a multivariant phase show such behavior, which 
cannot be reproduced with the RFIM (see section |l|. In 
our model, this feature descends from the intrinsic dif- 
ference of the physical processes occurring in the two 
branches (transition from the state to the three vari- 
ant phase in the first branch, and the opposite process in 
the second branch). In order to study this feature more 
quantitatively, wc define an asymmetry factor A as: 



A- 



{dM/dH)i - {dM/dH)2 
{dM/dH)i + {dM/dH)2 



(4) 



where {dM/dH)i and {dM/dH)2 are the derivatives of 
the hysteresis cycle at the coercive fields (defined as the 
two fields at a height m = 0.5) in the two branches. 

As can be seen in Fig. [51 A is greater than zero for neg- 
ative values of A, while for A > the effect of the 'dipolar' 
term is screened by the disorder and the Potts term and 
A is essentially 0, irrespective of the value of a. We have 
cut the curve with ct = 3 in Fig. [8] at the value A = — 7. 
In fact, for the considered range of parameter, lower A 
hysteresis curves begin to show the plateaux explained 
above, and thus our definition of asymmetry loses sense. 

The effect of disorder on hysteresis loop properties is 
illustrated in Fig. [31 where we show the width W of the 
loops at m = 0.5 (coercivity) as a function of increasing 
correlation p, for different values of a. As mentioned in 
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FIG. 6: (Color on line) Hysteresis loop for various values of 
the system size L = 16, 20, 30, 40, 60. The model parameters 
are: A = -25, a = 3, p = and AH = 0.05 
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FIG. 8: (Color on line) Hysteresis loop asymmetry as a func- 
tion of the parameters A, for ct = 3, 4, 5, 6. The system pa- 
rameters are: : L = 16, AH = 0.05 and p = 0. Each point 
represents an average over 800 realizations of the disorder. 
Solid lines are a guide to the eye. 
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FIG. 7: (Color on line) Hysteresis loop area as a function of 
the parameters A, for cr = 3, 4, 5, 6. The system parameters 
are: L = 16, AH = 0.05 and p = 0. Each point represents an 
average over 800 realizations of the disorder. Solid lines are 
a guide to the eye. Error bars are not visible on the scale of 
the picture. 



FIG. 9: (Color on line) Hysteresis loop width as a function of 
the parameters p, for cr = 2, 4, 6. The system parameters are: 
L = 16, AH — 0.05 and A — —3. Each point represents an 
average over 800 realizations of the disorder. Solid lines are a 
guide to the eye. 



B. Microstructures 



the qualitative description above, both p and a decrease 
the width of the loops. 



As we have already mentioned (see section |TT]) , when 
the dipolar term is large enough compared to the Potts 
term, the transformed domains lose their spherical sym- 
metry and start to show a non-trivial microstructure. 
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The microstructure of a system is defined as tfie arrange- 
ment of the variants of the product phase. 




□ = 



(0,0,0) 
(1,0,0) 
(0,1,0) 
(0,0,1) 



FIG. 10: (Color on line) Saturation configuration for system 
parameters: L = 32, AH = 0.05, A = -72, cr = 20 and 
p = 20. Different colors correspond to different spin variants 
(see the legend). 

In Fig. [10] we can see an example of these three- 
dimensional microstructures. We represent the views of 
the yz, xz and xy surfaces, when the sample has reached 
saturation (fuUy transformed state). In this case (A < 0), 
the domains tend to be prolate. For instance green do- 
mains (corresponding to Si = (0,0,1)) tend to be elon- 
gated along the z direction both in the xz and zy planes. 
This effect can be quantitatively measured as will be ex- 
plained below. 

For A > 0, we observe the formation of oblate domains, 
as shown in Fig. [Til developing in the plane perpendicu- 
lar to the spin direction. This effect generates a sort of 
"chessboard" correlation as can be seen, for instance in 
the yz plane by observing the red and green domains. 



In order to quantify the shape of the domains in such 
microstructures, we have calculated the average linear 
size {Dx), (Dy) and (Dz), of the domains of the three 
variants x, y and z, along the three spatial directions x, 
y and z, at the saturation configuration. For instance, the 
average size matrix corresponding to Fig. [TUland Fig. [TTJ 
are shown in Tables [ij and [III In the first case (corre- 
sponding to the prolate domains) , the diagonal elements 
of the matrix are sensibly larger than the others, con- 
firming the growth tendency of domains along the orien- 
tation of each variant. As is quite obvious, for symmetry 
reasons, the diagonal elements can be averaged giving 
what we will call the average linear size in the parallel 
direction {D\^) and the off-diagonal elements can also be 
averaged giving the average linear size in the perpendic- 
ular direction {D±). For the case of Fig. [TOj (A < 0) we 
get (Dii) = 4.44 ±0.20 and {Di_) = 2.19 ± 0.05. For 
the case of Fig. [TT] (A > 0), the tendency of the sys- 
tem to form oblate domains is confirmed by the values 
(Dll) = 1.175 ± 0.005 and {D±) = 1.89 ± 0.03. 







{Dy) 




X 


4.35 ± 0.14 


2.21 ± 0.02 


2.20 ± 0.02 


y 


2.04 ± 0.02 


4.22 ± 0.13 


2.15 ± 0.02 


z 


2.27 ± 0.02 


2.28 ± 0.02 


4.75 ± 0.07 



TABLE I: Average size matrix for the saturation configuration 
of Fig. M\ 







{Dy) 


(D.) 


X 


1.172 ± 0.003 


1.927 ± 0.014 


1.903 ± 0.014 


y 


1.862 ± 0.013 


1.173 ± 0.003 


1.868 ± 0.013 


z 


1.876 ± 0.013 


1.887 ± 0.013 


1.181 ± 0.003 



TABLE II: Average size matrix for the saturation configura- 
tion of Fig. [TT] 
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FIG. 11: (Color on line) Saturation configuration for system 
parameters: L = 32, AH = 0.05, A = 3, ct = 1 and p — 1. 
Difi'erent colors correspond to difi'erent spin variants (see the 
legend) . 



For the sake of completeness, as a final microstruc- 
ture example, in Fig. [12) we show a system configuration 
with high disorder and A = 0. As expected, no domain 
asymmetry arises and every spin just aligns with its lo- 
cal random field. The values of {D\^) = 1.51 ± 0.14 and 
{D±) = 1.50 ± 0.20 are equal to within statistical errors. 





(D.) 


{Dy) 


(D.) 


X 


1.509 ± 0.082 


1.491 ± 0.080 


1.492 ± 0.079 


y 


1.517 ± 0.083 


1.511 ± 0.081 


1.504 ± 0.080 


z 


1.512 ± 0.082 


1.512 ± 0.082 


1.515 ± 0.081 



TABLE III: Average size matrix for the saturation configura- 
tion of Fig. [12] 



The formation of microstructures such as those in 
Fig. [TO] and [11] is quite clearly affected by the dynamics 
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FIG. 12: (Color on line) Saturation configuration for system 
parameters: L = 32, AH = 0.05, A = 0, cr = 40 and p = 40. 
Different colors correspond to different spin variants (see the 
legend) . 



due to the effect of kinetic constraints. In fact, when a 
domain of one variant starts to grow, it necessarily blocks 
the growth of neighboring domains of other variants and 
vice versa. Thus the first variant to locally break sym- 
metry will facilitate the nucleation of large domains, thus 
restricting the dimensions of the other variants. This ef- 
fect could be seen, for instance, by analyzing the decreas- 
ing field branch in Fig. [13) the formation of the domains 
of type X and y that cross the system block the growth 
of domains of the z variant. 

Moreover, with the help of the microstructure repre- 
sentation, we can analyze in more detail the bump for- 
mation due to finite size effects, discussed in section HVAI 
(see Fig. [T5)) : in fact, if the system size is finite there is 
the formation of domains spanning a whole system side, 
such as the y domains in microstructures labeled with 
4 and 5 in Fig. [T21 As already pointed out, these kind 
of slab domains are actually infinite due to the periodic 
boundary conditions and are thus very stable. As we 
can see in the figure, they keep existing up to high driv- 
ing field, giving rise to the existence of the plateau and 
disappear when the field overcomes a certain threshold. 

As we already mentioned in section [TTl the truncation 
of the dipolar term does not allow elastic effects to be re- 
produced, which would lead to more realistic microstruc- 
ture. In particular, it would be interesting to control 
the tendency for the different variants to exhibit a pre- 
ferred habit plane, as would be the case of a real cubic to 
tetragonal transition. We expect that including a next- 
nearest-neighbor dipolar interaction will allow for such 
an interesting property to occur. The model presented 
here is, therefore, a promising starting point for modeling 
materials with a phase transition from a single variant to 
a multivariant phase. 



C. Avalanches 

Another phenomenon that may be analyzed with our 
model is the avalanche dynamics. The hysteresis loops 
in athermal first-order phase transitions consist of a se- 
quence of jumps between metastable states. Such discon- 
tinuities are, in general, microscopic. However, for cer- 
tain values of the disorder one or more may become com- 
parable to the system size and then correspond to the ob- 
served macroscopic discontinuities in the ferromagnetic 
loops. In the magnetic case, microscopic avalanches can 
be detected experimentally by appropriate coils. They 
correspond to the so-called Barkhausen noise [l3|. In 
structural transitions, avalanches can also be detected 
typically by acoustic emission techniques [13]. Knowl- 
edge of the distribution of the number of avalanches along 
the transition, as well as their size and duration, is an 
important piece of information in order to characterize 
athermal FOPT. 

Good discrimination of individual avalanches in the 
simulations can only be performed in the limit of AH — > 
0. This will require a true adiabatic simulation algorithm 
which is not available at p resent, as opposed to the case 
of the standard RFIM [14!|. In our case, after a small (but 
finite) change AH some spins can be updated. We will 
consider all of them as being part of a single avalanche. 
This is an approximation and, therefore, we should keep 
in mind the fact that we are slightly overestimating the 
size of the observed avalanches due to some unavoidable 
overlaps. In the experiments recording avalanches the 
same problem occurs [2l|, . 

With this definition we can study, for instance, the 
average number of avalanches (Nav) as a function of the 
external field H. As an example, in Fig. fnf a) we present 
the behavior of this number for the case A = —8, cr = 4, 
p = 0, compared with the behavior of the average hys- 
teresis loop in Fig. [T3] (b). (Nav) presents two peaks 
in correspondence with the two coercive fields and goes 
to zero (as expected) at the m — and m — 1 satura- 
tions. This kind of information is very relevant for the 
understanding of the measurements of acoustic emission 
in structural transitions using the pulse-counting tech- 
nique [23l |. 

More interesting information can be obtained by mea- 
suring the avalanche sizes S and computing the inte- 
grated probability distribution P{S), by analyzing all 
the avalanches in a single branch of the loop (the two 
branches must be analyzed separately since they are not 
symmetric). As a naive approximation, in our case one 
can define in our case the size S of the avalanche as 
the order parameter variation AM associated with an 
avalanche (i.e. when the field is varied by AH.) Never- 
theless, this definition imported from the standard RFIM 
should be carefully adapted to our multivariant FOPT. 
Inside an avalanche, in fact, one can distinguish between 
different kinds of processes taking place, depending on 
their effect on the order parameter variation. Let us focus 
on the decreasing field branch starting from the m — 
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FIG. 13: (Color on line) Parameters: L — 32, A — —12, a — 1.5, AH — 0.05 and p — 0. The configuration snapshots are taken 
for (1) m = 0.03; (2) m = 0.98; (3) m = 1; (4) m = 0.78; (5) m = 0.28 and (6) m = 0.18. 



phase up to the m = 1 saturated configuration: there are 
several microscopic possibiHties for a spin jump. A spin 
could jump from the state to one of the three variants x, 
y and z thus giving rise to a positive contribution to the 
magnetization change AM; it could jump from the x, y 
or z states to 0, causing a negative contribution AM < 0; 
or finally there is a possibility of a change from one vari- 
ant to another without contributing to the change in the 
system magnetization. These three possible updates will 
be called +, — and 0. Instead of only measuring the total 
size AAf of an avalanche, for each of them we will mea- 
sure the three quantities and which are the 
number of spin updates of each kind. Moreover, since as 
we have seen in section UVAi the hysteresis loops are not 
symmetric, we should separately analyze the avalanches 
during the decreasing field branch and during the in- 
creasing field branch. This gives, therefore, 6 possible 
distributions: P^°'""(n), P^°'"''{n) for the de- 
creasing field branch and P^^{n), P^^{n), PQ^{n) for the 
increasing field branch. It should be mentioned that in 
the increasing field branch the total number of events of 
the + and kind are much smaller than the number of 



— events, typically by 2-3 orders of magnitude. For the 
decreasing branch the — events are rare, but events 
are frequent, since an important number of transitions 
within variants may occur during the avalanches in the 
later stages of the transition. Of course such an optimiza- 
tion between variants cannot occur in the reverse process 
during the increasing field branch. 

In Figs.ll5land ll6l we show some examples of the prob- 
ability distributions for varying values of the two disorder 
parameters on log-log plots, respectively a and p. Actu- 
ally, we show only the P^°'""(n) and the P^^{n) distri- 
butions which account for the majority of the events. In 
Fig.fTHb) it is possible to notice another finite size effect: 
the slab domains of size corresponding to multiples of the 
system size L present the tendency to disappear abruptly 
(and thus the P^pin) presents some peaks in correspon- 
dence of values multiple of L) , as we have already pointed 
out in section HVBI 

An interesting result for the case of p = (see 
Fig. [IHa)) is that the distribution P^°'""(n) shows a 
exponentially damped behavior but seems to become a 
power-law for a certain critical value of a (which will 
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FIG. 14: (a) Number of avalanches Nav (arbitrary units) as 
a function of the external field. Nav is computed in each 
Ai/ interval over 600 loops, and (b): The related averaged 
hysteresis cycle. Simulation parameters: L — 16, AH = 0.02, 
A = -8, cr = 4, p = 0. 



be located around a ~ 3). Such a tendency has been 
well studied in the standard RFIM. The fitted value of 
the power-law exponent is 1.1, clearly smaller than the 
expected universal value 2.03 ± 0.03 [l^l, but this fea- 
ture could be due to finite size effect. Only a detailed 
finite-size scaling analysis and the study of the geometry 
of the avalanches [2^ (out of the scope of this paper) 
will reveal if the observed power law conforms to univer- 
sal behavior or not. Interestingly it seems that for the 
increasing field branch (see Fig. fW b)). the distribution 
is always exponentially damped, at least for the studied 
range of values of a. Therefore, the critical point for the 
increasing field branch, if it exists, would be located at a 
different (smaller) value of a. 



V. SUMMARY 

The analysis of microstructure formation in ferroic sys- 
tems undergoing a first-order phase transitions is an in- 
teresting issue both from a purely theoretical and an ap- 
plicative point of view. Microstructures arise since the 
product phase, arising from the balance of many ener- 
getic terms, may show energetically equivalent variants. 
Despite the interest in this issue, the models that have 
been used up to now for the study of the interplay be- 
tween disorder and athermal evolution (for example, the 
Random Field Ising model [Tlj) are not suitable for the 
analysis of microstructure formation, due to the equiva- 
lence of the variants of the product phase. 

In the present work, we have introduced a modifica- 
tion of the Random Field Potts model, which consists of 
adding a dipolar term truncated to the nearest-neighbor 
approximation, which represents a promising step to- 
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FIG. 15: (Color on line) Avalanche size probability distribu- 
tions P+°™"(n) (a) and P"*'(n) (b) corresponding to the two 
branches of the hysteresis loop, averaged over 600 disorder re- 
alizations and for four values of a as indicated by the legend. 
Simulation parameters are : L = 16, AH — 0.05, A — —8 and 
p = 0. 



wards the analysis of athermal transitions from a degen- 
erate to a multivariant phase. 

In our simulations we have chosen extremal updating 
in order to preserve the independence of the trajectory 
from the applied field rate. We have studied the depen- 
dence of the hysteresis loop shape on the hamiltonian 
parameters values. From a quantitative point of view, 
this has been performed by measuring the the loop area, 
its asymmetry and width. Our loops display a large area 
and asymmetry regime for very negative values of the 
'dipolar' term parameter A, associated with the forma- 
tion of microstructures with prolate domains, oriented 
along three equivalent spatial directions. On the other 
hand, for A > domains are planar (oblate) and the 
loops show low, almost constant values of the area and 
the asymmetry. 

We have also addressed the study of the avalanches 
in the hysteresis loop, distinguishing between the two 
branches. For certain parameter values, the probabil- 
ity distribution of the size of the leading process in an 
avalanche displays power-law behavior, which is the typ- 
ical result for athermal phase transitions in disordered 
systems. 
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FIG. 16: (Color online) Avalanche size probability distribu- 
tions Pf°^"(n) and P"*'(n) corresponding to the two branches 
of the hysteresis loop, averaged over 600 disorder realizations 
and for three valued of p. Simulation parameters are: L — 16, 
AH = 0.05, A = -8 and a = 8. 
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